{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(100, 2)\n",
      "(100, 2)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH0lJREFUeJzt3XuQVNWdB/Dvj3HUQawaI5jg4OwQdX3GgI6vnYorrMZn\nJphN+YgmWcsU5RoTjAqiW7UmVrlAqAQfRdyw0YpZMWrUIoomaInExF01M4BRIW6hgjJoIOIYgVFm\nmN/+MT3Tj7ndfd/nnHu/nyqq5vZcus+l6V+f+zu/c46oKoiIyB1jTDeAiIiCYeAmInIMAzcRkWMY\nuImIHMPATUTkGAZuIiLHMHATETmGgZuIyDEM3EREjtkriScdP368trW1JfHURESZ1N3d/VdVneDn\n3EQCd1tbG7q6upJ4aiKiTBKRTX7PZaqEiMgxDNxERI5h4CYicgwDNxGRYxi4iYgck0hVCRFRFMvW\n9GDhitexpbcPBzc3YfZZR2DG1BbTzbIGAzcRWWXZmh7c+Ogr6OvfAwDo6e3DjY++AgAM3gVMlRCR\nVRaueH0kaA/r69+DhSteN9Qi+zBwE5FVtvT2BXo8jxi4icgqBzc3BXo8jxi4icgqs886Ak2NDWWP\nNTU2YNqRE9AxfyUmz30CHfNXYtmaHkMtNI+Dk0RkleEByNKqkmlHTsAj3T3xDVju3Al8/DFw4IGx\ntTtNDNxEZJ0ZU1vKAnLH/JVVBywDB+677wa+9S3g0EOBDRviaG7qGLiJyHqxDVh2dgKPPz708003\nRWyVOQzcRGS9g5ub0OMRpH0PWO7eDeyzT/H4lVeAY4+NqXXp4+AkEVmv2oDl7LOOqP+Xn3mmPGjv\n2OF00AbY4yYiB3gNWPqaBn/xxcCDDxaPVRNsZXoYuInICZUDlnWJFH8+6CDgL3+Jv1GGMFVCRNmy\ne3d50L799kwFbYA9biKrcZW8gFavBk44oXj81ltABjcuZ4+byFLDq+T19PZBUZx0kucZgzVdd115\n0B4czGTQBtjjJrJGZe961+6B+CadZF1pagTIzCBkNexxE1nAq3f9wa5+z3O96plza3CwPGjfdFPm\ngzbAwE1kBa81qKsRgOkSYGi6ekNJbffLLwO33mquPSli4CayQJCp2wrkclOBZWt6RlYHvP+0C4HD\nDy/+sr8fOO44c41Lme8ct4g0AOgC0KOq5yfXJKL8qTalu5q8bSpQup3ZxgUV4ScHqZFKQXrcswCs\nT6ohRHlWbUp3c1Oj5/l521Rg4YrX0bd7oCxo9+47Dh3znjHYKnN89bhFZBKA8wDcCuDaRFtElEPV\npnQDKNs4FwiwRkeGNL71BjYumTlyPOv86/DrY6ZBErjzcKF23m+q5DYAcwDsn2BbiHKt1pRu2wNJ\noq69FquWLBo5POaah7Bzn7EA4r/zcGWH+bqBW0TOB7BVVbtF5PQa580EMBMAWltbY2sgUd4FXqMj\nQan3Rivqs9tuWD7ycxJ3HrV2mLflPQD85bg7AHSKyEYADwCYLiL3VZ6kqktUtV1V2ydMmBBzM4nI\ntNRnclYE7WWrN6OluQkCoKW5CfO+8rnYg6krO8zX7XGr6o0AbgSAQo/7elW9LOF2EZFlUuuNbt0K\nfPrTxeMFC4A5czADyacr/GzYYEMOnHXcRORLKr3Rb3+7PGhv2QLMmRPf89dRb8MGW9aPCbRWiaqu\nArAqkZYQkdUibx9WT8T1RuLoCdfbsMGWHDgXmSIiX2afdURypYkxBO24qkFqDQbbkgNnqoSIfJkx\ntQXzvvK5eAcIe3vLg/app4aaCVmrJxynancXaU+IYo+biHyLtTTx9tuBa64pHq9cCUybFuqp0uoJ\nJ3rXEQADNxGlrzI1Urk8a0CJ598LQm9aHDMGbiJKVwKbHqTZE7ZhQhQDN1FG2VBvXKa/H9h77+Lx\nkUcC6+NZt86WnnBaGLiJMsi6NTeefBI477zi8dNPA2ecEetL2NATTgsDN1EG2VJvDGBow95Nm4rH\n/f3AXgw9UfBfjyiDbKk3ztsmvmlh4CbKoLSqLKpSBcaMGf1YhpgcQ+AEHKIMqrfmRqJWry4P2vfe\nm8mgbXLNEva4iTLIWJXFOecAv/1t8bivD9h332Rf0wDTYwgM3EQZlXqVRY7y2abHEBi4iSg6S4J2\nWnln02MIzHETUXjr1pUH7csvNxq008o7e40hSOE1O+avTDzXzR43UUYl3vs84wzgmWeKx2+/DRxy\nSHzPH1CaeefSMYSe3j4IgOGvqzQmOzFwE2VAZZCeduQEPNLdk9zMSUtSI6XSzjsPjyF0zF85Km2S\n9EAlAzeR47ymty994W1UhtLYgomFQRuIJ+8c5i7FxEAlc9xEjvNKEVQLpZGCyTvvlAftE0+0JmgD\n0WvXw+bITWyuwMBN5LggwTh0MOnsBFpbi8d/+APw0kvhnishUXfoCbuLjonJTkyVEDmuWoqgdMAM\niBBMLE2NeIlSux425WFishMDN5Hjqm0i8M8ntODZP2+LFkwcCtpRRcmRpz3ZiYGbyHGJ9Ph27gTG\njSt/LMNBG7BnP0k/GLiJMiDWHt8ddwCzZhWPly4Fvva1eJ7bYi7tosPATURFlamRPXtGL8/quFol\nf67sosPATURDcpDPtm5Lt5Cy9VVKlGHL1vSgY/5KTJ77RLzrYezZk4ugDYQv+bMNe9xEFhu+rU9s\nPYzf/x447bTi8V13AVdeGf75LGd6Oda4MHATWaQ0/9o8thE7Ph5A/+BQuI59CvvJJ5dPotm1C2hK\naWszQ0wvxxoXpkqILFE55fqDXf0jQbua0D1FkfKgrZr5oA0Y3tItRuxxE1nCK/9aT6ieYk7y2V5c\nKvmrhYGbyBJBe8+Be4qrVgHTphWP77gD+M53qp6e1HreJndHB9wp+aulbuAWkX0BPAdgn8L5D6vq\nzUk3jChvquVfSw0PULYEDXiVvezt24EDDqh6elJlc1kpxzPNT477EwDTVfXzAKYAOFtETkm2WUT5\n45V/bWwQNDc1jqx2t+iiKdg4/zw8P3d6+KCtWjNoA8mVzWWlHM+0uj1uVVUAOwqHjYU/+UmKUaxM\n3ybbLJH8a5V8dr33IamyuayU45nmK8ctIg0AugEcBmCxqr6YaKsok7Jwm5z0F09s+deNG4HJk4vH\nn/0s8MYbAPy9D0mVzWWlHM80X+WAqrpHVacAmATgJBE5tvIcEZkpIl0i0rVt27a420kZEPdtcmIz\nCWu8Xlq7iEdy6aXlQXvNmpGgDfh7H5Iqm8tKOZ5pgapKVLVXRJ4FcDaAVyt+twTAEgBob29nKoVG\nifM22UTvPc1dxEPzUern531IqmwuK+V4pvmpKpkAoL8QtJsAnAlgQeIto8ypdps8RgTL1vQE+vCa\nCKLW52d91mf7TVckVTaXhXI80/ykSiYCeFZE/gTgjwCeVtXlyTaLssjrNhkA9qgGTjlUC5Y9vX2J\npS5MbArrS29vedBuaak5qca2dEXaKa8sqBu4VfVPqjpVVY9T1WNV9ZY0GkbZM7yZa0NlzxDBc921\ngmVSeWfbAh4A4JZbykv7fvc7YPPmmn8l6qa6cXJm3MAyoglMd21vb9eurq7Yn5eyYfLcJzzrSQXA\nW/PP8/UclTnuSgeMbcTYvffK3Ky/MpVfgIODox+zXMf8lZ5pm5bmJjw/d7qBFpkjIt2q2u7nXE55\np9TFURI2HCyveXCt5+8/2NWPD3b1A4h34NKa/GyK640k+WVl/biBpbg6IKUurpTDjKktaPEZ7DMz\nO29gIPWgnWQqw9pxA8sxcFPq4syxVhvw9FLZizM1KBb6dbu7gcbG4vHSpYmv7Jf0FHUrxw0cwFQJ\nGRFXysGrLnjnJwPo7esfdW5pL87ULM7Qr3vFFcA99xSPBwaABn9fWFEkncpgXXc4DNzkvMovAa+B\ny8penKnJNKFe1+D62WlMUbdm3MAhTJVQ5vhJxZgaFAv0uqrlQfvMM1Pf9ICpDDuxx02ZVK8XZ2qx\nI9+v+8YbwGGHFY+femoocKeMqQw7MXBTLs0+64i66RRjr9vZCTz+ePHY8Ca+UVIZJuveraq5jxkD\nNzkjzg+iqZ5k3dfN0H6QJpfxzcISwrVw5iQ5odqAo6mp2sNtSmPTA1eZnBXp4ozMIDMnOThJTrBt\ny6tYJ6Zs2lQetO+7ryxou7oIk8lZkVmfkcnATU6w7YMY2xfJ2WcDbW3F4+3bhzZCKHB5ESaTsyKz\nPiOTgZucYNsHMZYvEhFgxYrisccmvrbdaQRhspQw62WMDNzkBNs+iJG/SHzms2270wjC5PKxNi1d\nmwRWlZATbKsnDl1OuH07cOCBxeNvfhP4+c+rnm6i3jzu6h1T71GWZ2Syx00UQqge3e23lwft9etr\nBm0g/TsNl3PqecIeNznBxrrcQD26kKV+ad9pOLEhMjFwkxucDighgrapWX8u59TzhIGbAokjoIR5\nDicDyq5dwH77lT/mM2iburswtYYLBcMcN/kWR/4z7HPYVg5Y1+LF5UH7scd8p0dMlgDaVr1D3hi4\nybc4AkrY53AqoIgAV19dPB4YAL70Jd9/3eTdRdbL6LKCqRLyLY6AEvY5bCsHrCqG9UbSSldUS1ll\nuYwuKxi4ybc4AkqY56gMMIsummJfYNmzB9ir4uMUcpGoNJactbFKh/xjqoR8iyNdEfQ5nKgrXrq0\nPGjfeWeklf3SSFfYOpXe1QW10sYeN/kWR7oi6HNYXwbY0AAMDhaPP/oIGDcu8tMmna6wsUqHdwH+\nMXBTIHEElCDPYWOAGeHw+tk2lv1Z/yVtEaZKyGpBygBTvc12OGgDdlbpWP0lbRkGbrKa3wCTWi68\nu7s8aM+a5VzQBuws+3OuVt8gBm6ymt8Ak8pg2wUXAO0lO0u99x5w223xPX/O2XgXYCvmuB2Q5d2q\n/fCTE0/8Ntvx1EglGwcCnanVtwADt+Vs/IDZyO9gW6gvwYwFbcDegUBO/vGnbqpERA4RkWdFZJ2I\nvCYis9JoGA2xtd7WNn5uswPnwV97LRP5bC8cCHSbnx73AIDrVHW1iOwPoFtEnlbVdQm3jcAPmF9+\nbrMD9TInThzKYQ9btw446qjE2p82G8sByb+6gVtV3wXwbuHnj0RkPYAWAAzcKeAHzL96t9m+vwQz\nmBqplMa0ekpOoKoSEWkDMBXAi0k0hkbjSHt8fJWbJRC0bZzGbWM5IPnne3BSRMYBeATANar6N4/f\nzwQwEwBaW1tja2DecaQ9PjV7mZWb+AKxBW1bB5c5EOguUR//OUWkEcByACtU9cf1zm9vb9eurq4Y\nmkcUL8+qkp/9B/CTnxRP+vWvgc7OWF6vY/5Kz1RXS3MTnp87PZbXoGwQkW5Vba9/po8et4gIgLsB\nrPcTtIlsNqqXWZkaGRwc/VgEHFymJPjJcXcA+DqA6SKytvDn3ITbRZQ8r3x2jEEb4DRuSoafqpI/\nAIj3fzORSZ98Auy7b/ljMVeODKdkenr7IABKn52DyxQV1yohK6RWeXH11eVB+6abEgnawxN9gKGg\nPdzzYfUGxYFT3sm41CovKtMgfX2je94x8Jroo+CAJMWHPW4yLpVp/V757ASCNsABSUoeAzcZVy2g\n9fT2RU+ZeA04hkiNBEnlcECSksbATcbVCmiRNkP40Y+AMSX/xT/1qdBBO8jiVLVmu9o4i5Lcw8Ad\nET+I0XkFumGhUyYiwPXXF483bQLefz9U+4KmcqpNJwdg/4715AQOTkZg83Rmlwz/W13z4FrP3wfO\nDVekRjrmPYPZ7zdgRsiVGGrlrKut7+01nbxj/kor18Am97DHHUGYQTX20L3NmNqCljhywxVBu+2G\n5ZF7ttVev3lsY6AeNActKS4M3BEE/SCmtqGtoyKthLhqlWfQHhalSqVau1QR6IvbtkFLdiLcxcAd\nQdAPInezqS30UqMiwLRpI4dfvfSHZUF7WNiebbV2fdjX73l+tdexaYlediLcxhx3BEEXo+etcn2B\nlxr1KPV7d/5KIObNJ7zaNTyl3e/r2LREr617TpI/DNwRBP0gcjebmFWpz05jd5dla3qw85OBUY/X\nex1b1sBmJ8JtDNwRBfkgcruomLzxBnDYYeWPldRnJ92zrawmGnbA2Ebc/KVjMGNqS7jd5FPEToTb\nGLhTZNOtsrMqN/F99FHgggtGnZZkz9YrzQAAY/feayRo214myk6E2xi4U2bLrbKTLNnEt16awYX8\nMTsRbmPgJjdYErSB+mkGV/LH7ES4i+WAFKvYa4N37rQqaAP1y/psq9em7GHgtkDQYGfrxInYa4Nv\nuQUYN654vGiR8aAN1K83t6lem7LJ1y7vQXGXd/+8KhSaGhuqTjwJen7l300ypxnrjuaVvez+fmCv\n9DJ7Uf+tbK8qIfvEuss7JSvoQFbYga80Kh1iy+0aTo3E8W/F/DEliakSw4IGu7DBMY3p9pFzu4OD\nxoM2YNfSBLamxcgsBm7Dgga7sMExjUqHKLnd/73zF0BD8e++fvnV1pb7pYXriVA1DNyGBQ12YYNj\nGpUOURaJOvW73xw5POp7D2NGy3nGApQtVSE29fzJLsxxGxZ0IkTYiRNpzZSLukjUyKp+Bies2DKr\n0JaeP9mHgdsCQYNdmIEvK2fK1Vg/GzAXoNJY68TPc3M9EaqGgTtHrKl0WLsWmDp15PBXp16A2add\nMeo0kwEqqX+rIBUrcfX8WZqYPcxxU7quuKIsaOPdd9G4+M7cTFgJkrcOPWZQggOc2cQeN6WnSqnf\njM8MHeahVxg0bx215+/CglcUHAM3paNOfbY1aZyEpZ235gBnNjFVQsnaurU8aF9trj7bBmmvY2JL\naSPFiz1uik3lINh/bX4KRy9eUDzhzTeByZPNNdACaVf3TDtyAu574W3Px8ldDNwU2bI1PfjB46/h\ng13FXc+fv/Gfyk/KcS+7UpppoWf/vC3Q4+SGuqkSEblHRLaKyKtpNIjcMly1UBq0Ny44v/wkBm1j\nmOPOJj857p8DODvhdpCjSqsWxn2yqyxor5p8AiZXTKqph4sqxYs57myqG7hV9TkA21NoCzlouOd2\n5QsP49XbLhx5fOYF/4Z/ufAHgQIEa47jx00dsok5bork4OamUfnstjmPAyKBA4SrNcc2z0y0cqkD\niiy2wC0iMwHMBIDW1ta4npYsNypoF1IjzU2N+H7nMYEChIv52DQ2qIgqLzXyeRJb4FbVJQCWAENb\nl8X1vC6zuScWWX8/sPfeZQ9NvmE5WiJcp4uLKrl6l0BuY6okIS70xEK76y7gqquKx/feC3zjG3ir\n5JQwX1q2LKcahIt3CeS+uoFbRH4J4HQA40VkM4CbVfXupBvmusz2xCqnrn/8MbDPPmUPhf3ScjEf\n6+JdArmvbuBW1UvSaEjWZLIn5nM/yChfWq7lY128SyD3ca2ShGSqflY10Ca+mfzSqiKOpVeJgmKO\nOyFp98QSGwj91a+AC4v12fjiF4EVK2r+lbylD1y7SyD3MXCH4CdIppmvTWwgtLKXvWULMHFi3b/G\n9AFRshi4AwoSJNPqiSUyEBogNVLJxUFGIpcwcAdkY7VI7DnlCEF7GNMHRMnh4GRANg68xTYQ+tJL\n5UF7v/24sh+RhRi4A7KxWiSWhYTGjAFOPrl4vHYtsGNHTC00gysNUlYxcAdk42prkUvSRMp71qrA\n5z+fSFvTwpUGKcuY4y5hW7VIEKFzyjHks21Ubywi0+vIUOaxx10QpIc2Y2oLnp87HYsumgIA+N6D\na927FX/vvcwGbaD2WAR74+Q6Bu6CWj00L05/+C+5pKwe+7tXLMSy1ZsNNih+tcYigr7XRLZh4C4I\nWi3i7IdfBHjggZHDtjmP47HxR7nzpeNTrbEIGyuDiIJg4C4IWi3i5Ie/IjXSdsPykcec+NIJoNaA\nrY2VQURBcHCyIOg0bafW49ixA9h//7KH2jw28bX6SyeEagO2nJJPrmOPuyBoSZ2pssDAtcmzZpUH\n7YUL0THvGc9TrfzSSQBX9CPXiSZQSdDe3q5dXV2xP69t0i4pq1wnBRj6sqgadCqrRnbvBhobgz+P\nRVjGR1klIt2q2u7nXKZKIkh7PY5A66TUKPWztRa9nkxvB0cUAAO3Q3wNiA4OAg3lKRyv+mwXF4Gy\ncYEvIhNyl+N2ef2KutUQDz1UHrRnz87NpBqiPMlV4HZ60gzqDIj+4z8CF11U/MVHHwE//GHKLUwW\ny/iIhuQqcDs7aaagajXE8ZOA554rnqgKjBtnrJ1JsXGBLyITcpXjzsKt9qjcdIbXG6nk6qAqUdxy\nFbidmjRTz4YNwOGHF48XLgSuv95ce1Li4qAqUdysCdxp1OdmZsbclVcCP/1p8XjbNmD8eCNNYV01\nUfqsmIDjNSGksUGw39574cO+/lgDgvOBxqLUiMsTeYhs49wEHK9Bw/49it6+fgDxTrRw+lbboqAN\nsK6ayBQrqkr8DA66VP0Ru23byoP2VVcZD9pANgZ7iVxkReD2OziYy4CwYAFw0EHF4w0bgMWLzbWn\nBOuqicywInB71ed6yV1AEAHmzi0eqwKHHmquPRWyVFft8oxayh8rctyV9bnNYxux4+MB9A8W0wGu\nBoTQLMtne8lKXTUXryLXWFFV4sX56o+wdu0C9tuveHzmmcBTT5lrTw50zF/pWd/f0tyE5+dON9Ai\nyiPnqkq8OF39Edb99wOXXlo8fukl4MQTzbUnJzjISq7xFbhF5GwAtwNoAPAzVZ2faKvyqLERGBgo\nHg8Ojk6XUCIyNaOWcqHu4KSINABYDOAcAEcDuEREjk66YbkiUh60VRm0U5SlQVbKBz897pMAbFDV\nNwFARB4A8GUA65JsWJKsyZ9XbnrQ0gJs3px+O3IuK4OslB9+AncLgHdKjjcDOLnyJBGZCWAmALS2\ntsbSuCRYU0Hw8svAlCnF4xdfBE46Kb3XpzK5HFMhZ8VWx62qS1S1XVXbJ0yYENfTxs6KNbmnTCkP\n2v39DNpE5JufHncPgENKjicVHnOS8QoCB+qzichufgL3HwEcLiKTMRSwLwbwtURblSBjFQSqwJiS\nG5zJk4E330z2NX2yJudPRL7UTZWo6gCAqwGsALAewEOq+lrSDUuKkQqCjRvLg/YTT1gVtF3eh5Mo\nj3zluFX1SVX9e1U9VFVvTbpRSaq6b2NSPcx584Z618N27ADOPTeZ1wrBipw/EQVi7czJJKVWQeBA\nPtt4zp+IAstl4B6WaG63NGhPnAhs2RLP88aMswaJ3GPFsq4mJJbb7e0tD9qPPGJt0AY4a5DIRbnt\ncUfddsuzt755NdDZWTzpgw+A5ua4mx4rzhokck9uA3eU3K7X7MvxX+0E3lxdPMnCfHY1Seb8WWpI\nFL/cBu4oud3K3vrGBecXf3nUUcA6Z5dxiZXXF9zsh1/G9x97DR/29TOQE4WU2xx3lNxuaa/8gfuL\nW4vNPmcWg3YJr3RU/x5Fb18/a8aJIsht4I5Sz13aK3/xkGMBACd++7/xP6d1VvsrueQn7cSacaLg\ncpsqAcLndmefdcRICmDRFy7Doi9chqbGBsxjJUaZaumoSqwZJwomtz3uKFKffekor3SUF9aMEwWT\n6x53FFy/ub7KUsPmsY3Y8fEA+geLFTesGScKjoGbElX5BcfyQKLoGLgpVbxTIYqOOW4iIscwcBMR\nOYaBm4jIMQzcRESOYeAmInKMaAKr2InINgCbYn9iYDyAvybwvCZl7Zqydj1A9q4pa9cDZOOa/k5V\nJ/g5MZHAnRQR6VLVdtPtiFPWrilr1wNk75qydj1ANq+pFqZKiIgcw8BNROQY1wL3EtMNSEDWrilr\n1wNk75qydj1ANq+pKqdy3ERE5F6Pm4go96wL3CJyj4hsFZFXq/z+dBH5UETWFv78e9ptDEpEDhGR\nZ0VknYi8JiKzPM4REblDRDaIyJ9E5HgTbfXD5/U48z6JyL4i8pKIvFy4nh94nOPM+wP4viZn3qNh\nItIgImtEZLnH75x6jyJRVav+ADgNwPEAXq3y+9MBLDfdzoDXNBHA8YWf9wfwfwCOrjjnXAC/ASAA\nTgHwoul2R7weZ96nwr/5uMLPjQBeBHCKq+9PgGty5j0qafO1AO73ardr71GUP9b1uFX1OQDbTbcj\nTqr6rqquLvz8EYD1ACrXNv0ygF/okBcANIvIxJSb6ovP63FG4d98R+GwsfCncvDHmfcH8H1NThGR\nSQDOA/CzKqc49R5FYV3g9ukfCrdCvxGRY0w3JggRaQMwFUM9oFItAN4pOd4MB4JhjesBHHqfCrfg\nawFsBfC0qjr//vi4JsCh9wjAbQDmABis8nvn3qOwXAzcqwG0qupxAO4EsMxwe3wTkXEAHgFwjar+\nzXR7oqpzPU69T6q6R1WnAJgE4CQROdZ0m6LycU3OvEcicj6ArarabbotNnAucKvq34ZvAVX1SQCN\nIjLecLPqEpFGDAW5par6qMcpPQAOKTmeVHjMSvWux9X3SVV7ATwL4OyKXzn1/pSqdk2OvUcdADpF\nZCOABwBMF5H7Ks5x9j0KyrnALSKfEREp/HwShq7hfbOtqq3Q3rsBrFfVH1c57TEA3yiMjJ8C4ENV\nfTe1Rgbg53pcep9EZIKINBd+bgJwJoA/V5zmzPsD+Lsml94jVb1RVSepahuAiwGsVNXLKk5z6j2K\nwro9J0Xklxga7R4vIpsB3IyhgRWo6n8C+CqAfxWRAQB9AC7WwpCyxToAfB3AK4WcIwDcBKAVGLmu\nJzE0Kr4BwC4Alxtop19+rsel92kigHtFpAFDweshVV0uIlcCTr4/gL9rcuk98uT4exQaZ04SETnG\nuVQJEVHeMXATETmGgZuIyDEM3EREjmHgJiJyDAM3EZFjGLiJiBzDwE1E5Jj/Bx4o3uDcnRmmAAAA\nAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x290d2534358>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def pca(data, topN=5):\n",
    "    data_mean = data - np.mean(data, axis=0)\n",
    "    data_cov = np.cov(data_mean,rowvar=0)\n",
    "#     data_cov = data_mean.T.dot(data_mean)/(N-1)\n",
    "    data_vals, data_vects = np.linalg.eig(data_cov)\n",
    "    data_index = np.argsort(data_vals)[:-(topN + 1):-1]\n",
    "    data_vects = data_vects[:, data_index]\n",
    "    data_recon = data_mean.dot(data_vects).dot(data_vects.T) + np.mean(data, axis=0)\n",
    "    return data_recon\n",
    "\n",
    "N = 100\n",
    "dataX = np.linspace(2, 4, N)\n",
    "dataY = dataX * 2 - 4\n",
    "dataX = dataX + (np.random.rand(N)-0.5)*1.5\n",
    "dataY = dataY + (np.random.rand(N)-0.5)*1.5\n",
    "data = np.array([dataX,dataY])\n",
    "a = pca(data.T, 1)\n",
    "print(data.T.shape)\n",
    "print(a.shape)\n",
    "plt.plot(dataX, dataY, 'o')\n",
    "plt.plot(a[:,0],a[:,1], 'r')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXt8VOWd/9/PTDIkgYAQKFCRgKgggZAgatCKFkUta3fr\nDbURQVRKaF39/X5qa6mulx9u23V369oKUhVYyEvxunZbf10FtV6IVYRgFa/lVjQCgqjILSTf3x8n\nZzKXc86cmcwtyff9ep1XZs55znOemcx8zjPf53sxIoKiKIrSeQjkegCKoihKcqhwK4qidDJUuBVF\nUToZKtyKoiidDBVuRVGUToYKt6IoSidDhVtRFKWTocKtKIrSyVDhVhRF6WQUZKLT/v37y7BhwzLR\ntaIoSpfkzTff/ExEBvhpmxHhHjZsGGvWrMlE14qiKF0SY8wWv23VVKIoitLJUOFWFEXpZKhwK4qi\ndDIyYuNWFEXxorm5mW3btnHgwIFcDyXrFBUVMWTIEAoLC1PuQ4VbUZSss23bNkpLSxk2bBjGmFwP\nJ2uICLt27WLbtm0MHz485X7UVKIoijfV1TB3LjQ1pa3LAwcOUFZW1q1EG8AYQ1lZWYd/aahwK4ri\nTWMjPPggHH10WgW8u4m2TTpetwq3oiiJOXQIDhzIiIAryaPCrSjdhPp6GDYMAgHrb319Cp3YAn7/\n/XDppWkeYX5w6623snLlSsAKJvzss89yPKJ4dHFSUboB9fUwezbs22c937LFeg5QW5tER6EQBINw\n5ZVwyy1pH2c+cMcdd+R6CAnRGbeidAPmzWsXbZt9+6z9vgiFoLgYrr4aNm6E3/wGBg1K+zi9aGiA\nf/5n62862Lx5M8cffzzXXHMNFRUVnH322ezfv5+ZM2fy+OOPh9vde++9jB8/nrFjx/Lee+8B8PXX\nXzNr1ixOOukkqqurefrpp9MzKJ/ojFtRugFbtybeX19vCfnWrTB0KMyf3zYbr6qCU06xZtg+xdq1\nLweuv95a//Tiiy/grbegtdUy9VRWQp8+7u2rquBXv0o8zg8//JCHH36Y3/72t0ybNo0nnngirk3/\n/v1Zu3Yt9913H3fffTcPPPAA8+fPZ/LkyTz00EPs2bOHk046ibPOOouePXsmvmgaUOFWlG7A0KGW\necRpP1hrjQsXgoj1PMqUsm5dUtdKm1kmgi++sEQbrL9ffOEt3H4ZPnw4VVVVAJxwwgls3rw5rs0F\nF1wQPv7kk08C8Oyzz/K73/2Ou+++G7DcG7du3crxxx/f8UH5QIVbUboB8+dHiylY1o+9eyHWO+1r\nivkrI7h833LmzatKWmy9zDJOffmZGTc0wJlnWmujoZB1c5g4MblxOdGjR4/w42AwyP79+13bBINB\nDh8+DFiBNE888QQjR47s+CBSQG3citINqK2FRYugvNwS6rIya3a9a1d82xIOMIZ3aKSa/94yJrEd\nIwY/ZplkmTgRVq2CO++0/qZDtDvCOeecw7333ou0/URZl+Svko6iwq0o3YTaWti82TI19OoFzc3u\nbU3bNoZ3rMjJMf4F3Da/+N3vl4kT4eabcy/aALfccgvNzc1UVlZSUVHBLVn2sDH2HSOdTJgwQbSQ\ngqLkL4FAuz07krVUU42LQJeWwpdfJuw71sYNUFJizfhtU8m7776bNXtwPuL0+o0xb4rIBD/n64xb\nUbohTrNfT9EGOPZYX33HmmXKy6NFW+k4KtyK0sXwEyE5f741C47EVbQDAZgxA/7wB99jiDTLbN6s\nop1uVLgVJUekJQTdoc/Zsy0XPJF2V7zYvp1mxa60tsKmTVkPuFHcUeFWlBzgV2CTJZkIydhZsSuD\nB8OKFRm50SipocKtKDmgwyHoLqTVFS8YhJkzYe1a6lcNysiNRkkNFW5FyQGZ8HWGNLni2YK9bRss\nXgyDBmXsRqOkhgq3ouSATPk6Oy06lpRY+x2JrG5TVWU9jhBsm0zdaHLJKaec4rg/NslUMtx2223h\nMPhMosKtKDnASWALC60Q9I7YkJN2xWtstHJrDx9uRbb87GeOi5CZutHkktWrV+d6CCmjwq0oOcAp\nBN0YKwS9ozbkpF3xWlvh4MF2AXeobpP0TL4T0KtXL8DKO/KjH/2IkSNHctZZZ7Fjx45wmzfffJPT\nTz+dE044gXPOOYemtvflt7/9LSeeeCLjxo3jwgsvZF+sHSnDaJIpRckRtbXtojpsWHzeEK/ETB3F\nTru6OXKnLeALFsDatfDaa1FjBf+pWpPh+j9eT+OnyeVDSUTVoCp+da6P7FXAU089xfvvv8+GDRvY\nvn07o0ePZtasWTQ3N3Pttdfy9NNPM2DAAFasWMG8efN46KGHuOCCC7jmmmsA+NnPfsaDDz7Itdde\nm9bX4IUKt6LkAdm0IdfXw5jp1dwkHkk/HAraRt5ouhIvvfQSl112GcFgkG9+85tMnjwZgPfff5+3\n336bKVOmANDS0sLgwYMBePvtt/nZz37Gnj172Lt3L+ecc05Wx6zCrSh5QKJ82elk3jzYLI2MZEPU\n/hYCBItCMGtWVsuS+Z0ZZxsRoaKiggaHkjszZ87kv/7rvxg3bhxLlizhxRdfzOrY1MatKFnCK4Al\nmzZkexZfxCEApG3bwPGWeSQHZclyyaRJk1ixYgUtLS00NTXxwgsvADBy5Eh27twZFu7m5mbeeecd\nAL766isGDx5Mc3Mz9TlwZlfhVpQskChSMp2JmRJFOP6loDrquZ3CtYJ34R//MYVX17k5//zzOfbY\nYxk9ejRXXHEFE9vyxoZCIR5//HF+/OMfM27cOKqqqsKeKHfeeScnn3wyp556KqNGjcr6mH2ndTXG\nBIE1wMcicp5XW03rqijRDBvmbAopL08Qbp4kflKqxtqvDxKihSB/O/NKRi73X1eyI2ha1+yldb0O\neDeJ9oqitJGtxcdkIxwFCNJC07emMXKZsw93Mmg+k+zgS7iNMUOAvwMeyOxwFKVrkq0AFs8bhB0l\nGYEBCmhhxOplcOmlHbq2kznoyivhtts61K3igN8Z96+Am4DWDI5FUbos2Vp89LxB2FGSEUgoBMXF\nMGcOrFjRoWs7zfabm+GOO6xiv0r6SCjcxpjzgB0i8maCdrONMWuMMWt27tyZtgEqSlcgU1VhYk0T\nU6c63yDW0rYg2WrNvWxPEnPssZaqpsGTxG22LwJZ9pbr8viZcZ8K/L0xZjPwCDDZGLM8tpGILBKR\nCSIyYcCAAWkepqJ0fpINRU9kL3YyTSxdahWrib1B9NsSHZloe5Lwzjvwgx/4Gn+i8bjN9o2BM87w\ndQnFLyLiewPOAH6fqN0JJ5wgiqKkzvLlIiUlIpYkW1tJibXfprw8+ri9lZc7dOjU0N7Gjk3LeJYv\nFzEmuk0gIPJP/xTf34YNG5J8R7oWTq8fWCM+tVj9uBUlD/HjHZI2TxUf/oh+xrN5c3TleGOsAMyu\nsjh56623snLlylwPA0gyAEdEXpQEPtyKonQcP6Lsy1Oluhr693e/0HHHwSuvRO1yMokkGs+TT1oZ\nYcFKTztnDrz6avpEOx/cDO+44w7OOuus7F/YAZ1xK0oeECtM/fo5t4sUZV+eKo2NUWkH7UVJjjsO\n1q+H99+HysqocThFeHqN5+WX4cIL2/e1tlr7J3rksEqGTNXn3Lx5M8cffzzXXHMNFRUVnH322ezf\nv5/GxkZqamqorKzk/PPP5/PPPweiCyz85Cc/YfTo0VRWVnLDDTcAsHPnTi688EJOPPFETjzxRF59\n9dWODdALvzaVZDa1cSuKf5zsx4WFIqGQt03ZPre83LItl5fHH4+1Z7e2GZ4/+tYVMuHIT+LOc7Ob\nl5U527j/6Z/ibdrFxSKrV3u/5mRs3EnZ8pNg06ZNEgwGZd26dSIicvHFF8uyZctk7Nix8uKLL4qI\nyC233CLXXXediIjMmDFDHnvsMfnss8/kuOOOk9bWVhER+fzzz0VE5LLLLpOXX35ZRES2bNkio0aN\ncr12R23cmh1QUXKMm/9zWRn06uWd/9oz1Wp1ddwuA9DaytGv/Cd3s5kz+FN4BgvuJpHdu2HZsvh8\n3D//eXubQADOOssyj6Rrtu01pnREnQ4fPpyqqioATjjhBP7617+yZ88eTj/9dABmzJjBxRdfHHVO\nnz59KCoq4qqrruK8887jvPMs6/HKlSvZsKE94+KXX37J3r17wwUb0okKt6LkGKccJmCJ5Wefpdhp\ncTEcOOB46DBBlnE5N9OuuvZCo1d62dibxG23wdtvW4+DQQiF0i/a9rUzlfK2R48e4cfBYJA9e/Yk\nPKegoIDXX3+dVatW8fjjj/PrX/+a559/ntbWVl577TWKioo6PrAEqI1bUXJIfb1jzQLA3a7sixjR\nDtu2R47kBN5kFkvYTnTAzdat/iM8H30Ubr/dehwKwTXXwKpV6RdtyG7K2z59+tC3b19efvllAJYt\nWxaefdvs3buXL774gqlTp/Lv//7vrF+/HoCzzz6be++9N9yusTG9VX0iUeFWlBwyb160C10kX32V\nPu+JcMDN+++zqIdz6lZ7Vp0owvPgQbjkkvbnLS3W+udll2XG6yNTUaduLF26lBtvvJHKykoaGxu5\n9dZbo45/9dVXnHfeeVRWVvKtb32Lf/u3fwPgP/7jP1izZg2VlZWMHj2ahQsXZmaAoIuTipJLYgNW\n0rYA59ZhQYE8/uumhME0btx/f/RCZDBoLaT26JFcfxqAowE4ipIT0uFbnMhOm/aak62tXPjDQSnN\nYBctio6Ov+EGuPNOy6Rz8GB0W69UskrHUeFWlBRIl2/x/PnuNm5IzwJc2BLTty8891zK/cQG0zQ2\nws03w44dzu0zUehYsVDhVpQUSLZgATjP0Gtr3W3cYAm775m9HSU5YwYQkQGwb19r5XD3bpg82fGm\nM326dQNx6//226GpKXqfHXSTrVzjSgR+bSrJbGrjVro6brZpY5zbOwXZGCNSV+cd9OInuVOYmECb\n/QHkv2c/EBeg43Y9t/737YtvM3q092tTG7c3HbVxq3ArSgokG83n1t4Wbzfh60gGwFaQZoKymCtk\nEJ+E+/USbaf+nY7PmhV96YQRnDGocOvipKJknWR9i72KDDzzjLu7m++oweLiuDZ2WbIZ/CePYJUl\n27fPCpZJhN3/t7/tfPyRR6JNKsnmGlc6hgq3oqRAsr7FXvberVvdhc+X/dglSlKAZoIsYSaX0F6W\nrKUl/qbj1P9997lXrlGvkdyiwq0oKZLMLNPLe8RL1H3N7F2iJL+gD+fwR2axOCpK0r7JlJdbz2PH\nVVIC554LP/yh83GbLVuylGbVLnIcuzqaRkSE1tZOVFLXr00lmU1t3IoST11d/KKmn8CXZDMARtq4\nv6RnwuvF9j99env7UEhk4MDENnG/ATw2Sdm47YEUFVlv4ief+D/Xg02bNslxxx0n06dPl9GjR8uS\nJUukpqZGqqur5aKLLpKvvvpKRERef/11mThxolRWVsqJJ54oX375pezfv19mzpwpY8aMkaqqKnn+\n+edFRGTx4sVy/vnnyznnnCPHHHOM3HjjjY7X1sVJRelEuIlwsot7UXgo6sEevZLqd/Xq6C4CAZGL\nL/a3qJlMlGfSwh15J0mTgG/atEmMMdLQ0CA7d+6U0047Tfbu3SsiIj//+c/l9ttvl4MHD8rw4cPl\n9ddfFxGRL774Qpqbm+Xuu++WK6+8UkRE3n33XTnqqKNk//79snjxYhk+fLjs2bNH9u/fL0OHDpWt\nW7f6ev3JCLeaShQlzXj5XTuZVzoUzOOQujXMiBGEXn/Vtzln0SKIyadEIAD/639F2/PdyErAzaFD\nlmno/vvh0ks73F15eTk1NTW89tprbNiwgVNPPZWqqiqWLl3Kli1beP/99xk8eDAnnngiAL1796ag\noIBXXnmFyy+/HIBRo0ZRXl7OBx98AMCZZ54ZTv06evRotrilf+wAKtyKkkZSEeFUgnnCRGSgs23b\njBhhVbf56KOo6jZe/PjHVjh7c3P7vmAQfvMbK+Nf5A3Hto3HkpWAm1DIWoydMwdWrEjcPgE9e/YE\nLMvDlClTaGxspLGxkQ0bNvDggw+m1GdsqtjDhw93eJyxqHArShpJRYRTLhQQM9sWwIRCcPbZMGBA\nwrHaNDTAL38Zv/+aa9oLLESSzTSrYWzBvvpq2LjRuqMMGpT4PJ/U1NTw6quv8tFHHwHw9ddf88EH\nHzBy5Eiampp44403ACsz4OHDhznttNOob7sbf/DBB2zdupWRI0embTyJUOFWlDSSSISdzCgphYxX\nV0fNtqHty3zoUNJmhMWL4/cFAnDFFc7ts51mlaqqjAm2zYABA1iyZAmXXXYZlZWVTJw4kffee49Q\nKMSKFSu49tprGTduHFOmTOHAgQPMnTuX1tZWxo4dyyWXXMKSJUuiZtoZx68xPJlNFyeV7opXpKNb\naLhX5KQrbiuEwaDIjBkiTU2+xnv//c4Lj7GRkelGIyd1cVJRsorX4qOXGcHNjBIZOQmWbXnfPitX\nlGPiJ68VwpYW2LTJ16z0nHMsu3bsmCA+MlLJM/wqfDKbzriVroqfhEpurn2JElM59e14DS+fvMGD\nfc22Tzklve59yaIzbp1xK0rW8LP4mGr4ulPfbteIROwHkybB2rUJZ9u33QarV3s2ATSfdj6jwq0o\nSZCyBwiJvTES9bF1K45+22HDyUsv+TKR/PznCZsAmk87n1HhVpQk6EjRgETeGIn6+EtBvCdJeLYd\nCsFjjyUcw8knx5cZcyLj7n1Kh1DhVpQk6KgPs1diKqe+I69R0Rwv2gZgyhSrs4su8rz25ZfD66+7\nHw8Gs+Tep3QYFW5FSYJUfZj9lB+L7Bva82YHg/DKPg8TycqVCf22Fy3y9hIpKYGlS5PPp52Ogsm5\nolevXh3uIxgMUlVVRUVFBePGjeNf//Vfs5Nl0O8qZjKbepUoSjuplPaKPc/T/WPgQE9PktWrRQoK\n4k8LhVJMatXB1yWSmlfJJ598IpMmTZImnz7qiejZs2da+9i+fbuceeaZcuuttyY8L+PZAYEi4HVg\nPfAOcHuic1S4FaWdZMucOZ3nKdzV1Z79XH11/Cm1tbl7XSKpCXddXZ0EAgGpq6tL+lwnbNFtbW2V\nG264QSoqKmTMmDHyyCOPiIjI3Llz5emnnxYRke9973vhbIAPPvig/PSnP43qw+avf/2r9OvXT1pb\nWz2vnQ13wIPAZBEZB1QB5xpjajIz/1eUrkeqnih+PFUEaFjXw9VMsWgRPPBA+/Ng0IqIX748cd+J\n6IiHTTIUFxdjjGHBggW0trayYMECjDEUO5RrS4Unn3ySxsZG1q9fz8qVK7nxxhtpamritNNO4+WX\nXwbg448/ZsOGDQC8/PLLTJo0ybGvo48+mpaWFnbs2JGWsbmRULjbbgZ7254Wtm3icYqiKBGk6omS\n6LiVDdBwPk85ZiFsaIC6uvbnxrgnjkqFjnjYJMPGjRv5/ve/T0nbym1JSQm1tbVs2rQpLf2/8sor\nXHbZZQSDQQYOHMjpp5/OG2+8ERbuDRs2MHr0aAYOHEhTUxMNDQ2ccsopabl2qvhanDTGBI0xjcAO\n4DkR+XNmh6UoXYdUPVHcvEzs9K2f04fJrAyXJYsM0mlogGuvtRYbbQoK3BNH2SSz2JitLIGDBw+m\nd+/eHDhwgKKiIg4cOEDv3r0ZlIFkU5EceeSR7Nmzhz/+8Y9MmjSJ0047jUcffZRevXpRWlrqeM7G\njRsJBoN84xvfyOjYklp0BI4AXgDGOBybDawB1gwdOtTbuKQonYQOVabpQD92eztvlGCVImsFkT59\n5AxWuYbPr14tUlgYva+gwEooleiayS42pvr+JGvjPv/882Xu3LnS2Ngoc+fOlfPPPz+p852w7dNP\nPPGEnH322XL48GHZsWOHDB06NLwAOmPGDDnqqKPkww8/lNdee02GDBki119/fVwfIiI7duyQKVOm\n5MfiZNwJcCtwg1cbXZxUugId8ZroyDXLyqxrtYLsJyQX8bAcxsieHiUiq1aJiPfC4KmnRu8bPdoS\n80QimyizYTpuYDb5kKsk0eKkiMgDDzwggwcPFhGRQ4cOSUlJiTzxxBPh44FAQMaNGyejR4+WyspK\n+Zd/+RdpaWlJeO1seJUMAI5oe1wMvAyc53WOCrfSFeiI10QqxN4oImfZ+wnJD/rVu7a1byrnnRc/\n3pNO8ncTckuCZbdN5w0sH4Q7l2TDq2Qw8IIx5i3gDSwb9+/TZqtRlDwlW14TNpFJpr7G8pgwbVsP\nDrFgdy306AGPPOIYCDR5Mvze4Zt51VX+kmO5LSraaWa9zlWyix+vkrdEpFpEKkVkjIjckY2BKUqu\nyYTXhNfin11T9muKKeFA1Hm2gHPoEFx2GRAdPv/ww/CHP8Rf76abLC8SPzcht8XGlpbE5yrZRUPe\nFcWFdHtNeBUSrq9vr48QK9pxxPgvNzTAjTdafUZy003wi19Yj91uNoFA+00EnMP5M1UcWGIH3E1I\ny+v2a1NJZlMbt9JViF2Uq6tLfZHOy2buJ0qy1XYRWb8+3Ofq1SI9ekQ3DQREbrop/nW4FWlIZLfO\nxCLtxo0bZefOnQkjDLsara2tsnPnTtm4cWPcMZKwcRd0XPoVpXuwdy88+KBlrYD2GTNYZou5c60Z\nakuLZReePRvuu6/9fD/mirXEJ5OKoroaKivDT3/xi+g0rSedBL/6FUycGH2anTRq3jzreoFAvAnE\ntlvHJpiKPXfoUOtXR0eyBw4ZMoRt27axc+fO1DvppBQVFTFkyJAO9WEkAz9XJkyYIGvWrEl7v4qS\nTWzThltVGpvycpg6FRYsiD9WV9cu3sOGtduxY88H2LAl3rZtI0BrIEjw423hYgkNDXDqqdEmkjlz\nnMcRSyAQb1oBy0SSjeR2SjzGmDdFZIKftmrjVhQXvEqJRbJ1qzXTdiJyv5fNfP78xLbtYNW4sGgv\nWgTf/W60+AaD7pGRsYui/fo5t9OqN50DNZUoigt+vSaczA42kfs9TQ4OJclsBDDf+EbYbWTRIqs6\neyTBoDWzjzWRgGXCWbiwXeS3bLEK5hQWQnNzezutetN50Bm3orjgd/bpJtrQXgwhITElyWLZ2fvo\n8Gz7kUfij59xBtx1V7ybYX19tGjbHDoEvXsnXxBCyRP8rmIms6lXidJZ8ArldvKmKCy0QtKNsXKI\neHlpgOWFYvdlh7I7emd4eJJsp79MGGLlznj5ZZFvfCO+aWRuEtv5xPaA8RpfJsP3leQgk7lK/Gwq\n3EpnwI+bm5ewe4WIB4PRou3livd2YZWnaA+kSYyxkkQFAolvFpHinahNpnOvKP5JRrjVq0Tptnh5\neWzenL7z3dqBc5SkjQBvMZYq3qJ3b/jyy8RjiiUY9DblOI1XyQ3qVaIoPuhoLhK/kZVe/SXyJBnG\nZiA10YbEog0aut4ZUeFWui0dzUXit+J7Ki52AuyjiG/xSjgUPlO4uQYq+YsKt9Itqa+3IiFjSdYl\nLjLR0+bNzl4ZbpVsysq8+36LKt4NVuJkzezRI/q5b+8VB/bsgf79/VW9UfIDFW6l22FHRO7aFb2/\nrCwzLnFOM/Ply+Gzr+OL3Urb9gKTqC1+ytHUcdNNVuh9ZH9Ll1p92vuSoaXFei9EcKxdqeQfujip\ndDs6uigZS319ink8PBS2xQQpkMNx+2tr/VVo91oQ9YMuWGYfXZxUFA/SWSDBK1WrJx6RkvsKSzlT\nno3b71e0wX3hNJF5xkYXLPMbFW6l25HOAgl+Kss44hIpKcCh5gB/YnLU/qoq/6Jt/wLYt6/d9m0v\nnN5zj7O9PRbNWZLfqHAr3Q632ejUqe7Vadxwm5lu2ZJ8XwK8x0hO46Wo/YFAdHpYLyJ/AYBlv7YX\nXG3zTWQdhl69rJwlkWjOkvxHk0wp3Q6nZE9Tp1oLfPbsOTbXthv9+sUvctrY4hnV19XFMGKEY/tW\nYAcDeZv2fNtFRfD8887Jo5xI9AsgNk1taytcfTU8+mj76yiOXzNV8gxdnFQUUl+w7N/fXbgd+9ri\nsCBZUMCBlkIelCu5k1vYzqDwocjyY37wyrM9dKjzaywrg/37owW9pESTTmUbXZxUlCRJdcFy927/\n19iwJXoqa+vrV6aUqfI7fsRvokT7+OMt0fYqMByLl/3e7bXs2qVV3DsbKtyKQuIFSzfxTGYRz6ly\nO0Cv5s95iguijg0YABs2JO+14hWGn+yCo3qW5C8q3IqCs+AVFlrRlcbA9OnO4ukWFRlLo/GuJbmJ\nYVHP5861/ibrteIVhp+si6B6luQxftMIJrNpWlelMxKZwrWsTCQU8k6JWl7efp5Xbu7ycvHsqBXk\nVWoc+/ZKzdrR12inqc1EFXcleUgiravOuBWljci8I716tVdzd8M2JdTWuhfYNcZ7cdMOcb+Ap6L2\nb9lizejdZr3GpC8s3W+yLCV/UOFWFAf82HcjRdXTRu7iXyfAYYJMZlXUoqTN7NmWm6JTZLxI8ouH\nXvZyP8mylPxBhVtRHPBj3506tX3RcsuWeIENB7IccM+5vZ5xcVGSNvv2WTNfN4/dZBcPU47yVPIO\nFW6lW+HXtc7PouOjj0ZHKYq0i3fY3HC396LkIUKex1ta3HNRpctLRL1HOh8q3Eq3wclUcOWVzrmo\nI+2+bjj5P4u0B+3U3l3tnZOEUJxt27GteMzmkyCdOVqU3KLCrXQbnEwFzc3uuahtu6+XeDsRnsE6\niLa9GPkhRzOBNxxt207YNwRjLPe94mLLRTHyZpPo14TfUmtKJyCR2wlwFPACsAF4B7gu0TnqDqjk\nI36qnke64tk4uct5uQqGz3dp0ALyIpN8jcXJ9dDJda+uzp9Ln1fVeiW3kIQ7oB/hHgyMb3tcCnwA\njPY6R4VbyUfKy/2JpDHx58YKXlmZ+7nLl4tIUZFjg1aQbQyWcQObUhJgt9fg5kceexNS8pdkhDuh\nqUREmkRkbdvjr4B3gSPTPvVXlAzjN8rRyeYb6y7nlqNEpM22HeNJYptIWghwG7cy945Bjr7T993n\n7VPttpDoVs1dFx67JkllBzTGDANeAsaIyJdu7TQ7oJKvRJYZ69cPvvgCDkdUCAuF4KGHEvsxe2YT\ndMoA2IZwTCQ2AAAa7klEQVT9bTNJfO/8XDcYdBZvLUHWechIdkBjTC/gCeB6J9E2xsw2xqwxxqzZ\nuXOn/9EqShaJnDnfc0+8t4ZfPXVb6FuLu/uf3fVL1z7me7x+rzt7ti48div82FOAQuB/gP/tp73a\nuJXOgJu92K9d2HGhz8N43gryBuM7nAfEbYFRFx47NyRh405oKjHGGGApsFtErvdzM1BTidIZ8Co6\n4JZ7xJNqb7/tT/kG1axnO4MyZsJIueK8knPSbSo5FZgOTDbGNLZtUzs0QkXJA9wCT0T814mMwkO0\nd1AWFm2IXzRMpliCGylXnFc6HX68Sl4RESMilSJS1bY9k43BKUom8fIySVr0qr1D2/9GeVSwTeRN\nI12Cq7lIug8aOal0OfzOXhOFte/bB9dd5/OiLrNtm+aInCSxi4bpElzNRdJ9UOFWuhT19Vb+kdh8\nJF7ivXmzeyKnXbsSz3wXLXLeL1iV21dwMee35SQJBqP9suvrnd37wN2cYgwUFFh//ZRR01wkXRC/\nq5jJbOpVouQKt4jGsjLv87yiKr28TFavFvka5yhJ25MkMrw9MirTKXzd7bpebe1QfqeqPVrJpvOA\nVsBRuiu7diW338bL39nL1HDBBfFFgMGebRtWMI1LWBHeHzn7dTKR2Pgxp4Sv1eYZYyfLKivTSjZd\nHRVupdvgZfKorU2+aO7ZZ8OnnzofsywvwmA+DS9Kxoqx1w0hVnD92qmbm62ya1rJpmujwq10KdzE\nFxIv9t1zj3P04dSp8YudDQ3w3HPwNc5lyQAkGOKGIStcZ79uN4Ty8uh29fXWtf2ii5HdAL82lWQ2\ntXEruWL5cnebsVPWP6fzI6MPnbL19ehhZeNzsm23tm0fBI+T/75rfcJrJUrFWlfnPx2tZgTs3KA2\nbqW7korJA9o9NqZPt54vW2aZGp55Jt62fPCgldAp1rZtB2F+yAgualnBJf+3MqF5xisTYH09LFzo\nnT8lHZVxlE6IX4VPZtMZt5JL/Mxk/bZ3m+2upcrTk+Rz+nR49uvl6WL/etD8JF0H0llIIZVNhVvJ\nNckImleyKbdjXqLdTEDGsD5l84w9Vi8TiZpDuh7JCLeaSpQuiR1Ys2yZ9Ty2PmMkXhGH8+fHmyOc\nFiSlbdtEOdWs420qgcTBL07h7tOnW9d0W5A0Rs0h3R0VbqXL4jcHiJu42sIZa2N28tu2tX03/cKi\n7cfe7OSfbV/PqTCCMTBnjrr5dXdUuJUui98cIG7JplpaYObM6H2HE3xljuUjwFponDHDupZXzhQ/\nrnvBYPvi5bJlVnkzpXtTkOsBKEqm8Jt0yZ69zpgRP8uNLGsGECRm+t2GAC0Y+qx/Balsn+3bNw57\nth95PbBm+265SmxaW1PMD650WXTGrXRZksm3XVubWBwTzbY3FI6HSstM0tHZfiSaJEqJRYVb6bIk\nm287kUB6zbYF+PDffx/el8xsPzK1rPplK35Q4Va6LH7ybUfOgP3Mft04HCrhwh86F0qIxGm/7QEj\nYtmw3QJyFMVGhVvp0iTKt711a3TUpJu5xC0niT0Hn3CoIcr84laNPdHsObIKvSaJUtxQ4Va6BW4z\n4H79ol0GD8R7+gHOLoA2eynhL1RGmV/cwtmh47UlFUWFW+kWuM2AwT3PtY3ToqRt1/6IEZxCQ3h/\npPkldvYMWsxXSQ8q3Eq3wG0GvHt34nOdFiVty8sAPgsH3Nhs2eIsxm6eJr7rWipKGyrcSrfByX6c\nyJOklXjjuC3jfwuWcxovOZ7nNJN28zTxU9dSUSJR4Vayht/q69lk7Fj3Y2updpDt9tn2kNZtbCyp\ndGjh7LPtdZNIVOQhH987JXeocCtZwW/ekGT6S4eQPfec+7FqGuP2hY0mPXsSWPmsa4V3iJ9hp1rX\nMt3vndIF8JtGMJlN07oqsXilTk2WZPNt2+fYqVPLytyrwUf2l7DUTAqvze26Xu9DOt87JX9B07oq\n+YbfSEI/XHedv3Bym9gZ665d7lXf7YXLvw2o9h7EY4+FHybjs+1W19JrNu6Wy0RrS3ZfVLiVrJBM\nJKEX9fXuousmZE7eHE7MmtW+cNlvS7yZBGBvvyGcO7iRwLSLwiaaRCXIIkmmLVj9uwUPaQ6Tbozf\nqXkym5pKlFhSMW844VXOy8104LfYbtR4XBq1YORFJiX1GjpSXszt9RqjZcq6GmjpMiUfSUd9RC8R\ndupv+XKrIrsf4QaRtwurrNLqLmXJtjFYBtLk29bsVKXd7w3Lq2J9hIld6SIkI9zGap9eJkyYIGvW\nrEl7v4oybJizzbesDD77LHpfbE5sPwgGgXg3wGCQxS3TuZl/ZjuDog4Z45zjpL7eyn/i9BUrL2+P\npnQi0dgTna90Powxb4rIBD9t1catdCrcFgKnTYt3D/Rr27ZZi7UgaYt2WG+POw7efJPbyxfHiTa4\n25rnzXMWbUi8sOg1dk31qiQUbmPMQ8aYHcaYt7MxIEXxorbWqlQTDFrPg0GYOBGWLo33c05UWcbu\nwybWbzs86/7gA5g1K+mMf17inGhh0etcTfWq+JlxLwHOzfA4FMUX9fWWSNslxlpa4Pnnnd0DY4U5\nktpaqx/bu2MdCdz/PvwwKY+Q+vrEVdq9gojchL28XEVbwd/iJDAMeNuv4VwXJ5VM4eVV4nebNcuh\nY5fGrSCHQiUi69f7HqOTB02kN0hdXWIvm3R54SidB9LtVeJHuIHZwBpgzdChQ7P1WpVuhl/XPhAp\nKIjfd/TRLh07CHYryFP8vUwY0pTUGN1uLsFgu/D6iYZMhxeO0nlIRrh9eZUYY4YBvxeRMX5m8epV\nomQKN6+SWIxxXhi8//72authiosdKygI8CbjOcm8mVSV9UDA+dqR3id+2ijdC/UqUfKeWPvu3Ln+\nkkbNnw+hUOL+nUTxvPMcRBvcy94Ax/Jh0hGKfqJE0xVJqnRPVLiVrOOU7W7BAn/Z72probQ0tev+\n5S8OOx1WEKVt+4ij+RavsHdvcpn4/HifpFqTUlGAxDZu4GGgCWgGtgFXJTpHFycVL/wuMHY0hN1p\nYTAOj0XJz+mT8sKgH/u02rCVSNDISSWfcbPvxhJr77WDavzYuJ2IizasrobGaN9te1ibKOcf+F1U\nWTKNVlQyidq4lbzGrx03sl2keSUVHM0QjfEZAO2gm6Fsi6slqWlUlXxBhVvJOk723VhsobUXMS+/\n3DtvRyylpQkCZarjA25s2/bn9OEsno07rguHSr5QkOsBKN0PW0DnzbNmsUOHwtSp8Mwz0c+vu849\n97aNMc6z8Lo6+MUvPE50mW0L0MRg/sTkqGO6cKjkE2rjVvKOVLL6AZx5piXiF1yQQLRd/LahfdYd\nRAgGLRv70KGWaGuouZJJkrFx64xbyTuSzepns3Klz4Yeot0KXIxVlqy11Xlx1O1Xgoq7ki1UuJW8\nI5VFwDlzOn5dAf5MDU9xEWB5vwQC7SK9dGn7DcX2Pbexfc9BxVvJPLo4qeQdGV0EdFiUBEu0D1PI\nBTwV3tfS0h4QtHBh4l8BXgWLFSWdqHAreYcfr5NYrrjCZ0OHRUmb9Yx1LJQA/vzOQV0Gleygwq3k\nHbF5r8vKoLDQvf2YMVYxhYS4uAA2E2AJM/guf0h5zDbqMqhkAxVuJS+prbWiFFtbrVqSP/2pe9uv\nvorfN3cuFBRYwl9QAH/rHx8lCZYLYADhaDaxnUGYuGKTEW09joG6DCrZQ4Vb6RS89JL7sVjzxNy5\n1sJhZJWco3Y5m0haA0Ge7HkFl7KC8nJvk8icOdFBPXV1/qrhKEq6Ua8SJe9paIAXXnA/HmueWLQo\n+nlrfM32MIHWFi4+YRMX/8mybffv7xz007Mn3Hef3xErSmbRGbeS1zQ0wCmneLeJNU/YM20bTwvH\nEUfAihUJx1FUlLCJomQNFW4lb2logO99z7tNWVm8eSKySPDaREWA9+yBQe2eJLt3Ozdz268ouUCF\nW8lLGhrg29+GHTvc25SUwD33xO+PrHJTjbv7HwCPPRb1VCvTKJ0BFW4lL1m8GA4ejN4XDMLAgYkX\nA++7z1o4/Jpi74scfzxcdFHUrqlTnZu67VeUXKCLk0re0dAAv/1t+/PCQrjqKivIxpe/NnDqqVCy\nwL2WJIMHw/PPx+1+5hnn5m77FSUXqHAreUfkYqQxlmhH5gXxw7x5EDsZtz39zNCh8LvfRdm2bdwi\nHzUiUskn1FSi5BUzZ7Y/DgYtb44rroivCp+oeK+T0Jq2ja1b4R//0fE8N1u2iL/rKko2UOFW8oaF\nC60MfADXXgt33gmrVsHGjfFV4d2qwNu8VVDluF/AckVxcQH0ypMSe91kbyaKkjb8VhVOZtMq70qy\nvPBCdKH11avbj3lVhXetju5Qtf0QRj781kyRpibPsdjV1xNds6Qken+yleAVJRKSqPKuM24lL/j2\nt9sfB4Pw4ovtz73sy06z7y17omuZCdBCkDUjr+CyTXcR+OYgzxmynSfFLTfJ1q3OxR40rauSLVS4\nlZwTKZDBIIRCcMYZ7fsS+VDHCuZdL98V3T9QQAsnvb+Muz++1Le5xcunWxcxlVyiwq3kjIYGqKho\nf/700+127Ui3Pz/5ubduBaqraSnrx6m3xiQrCYXYb4pZyBwuod22nWiG7HRdOwOgBuooOcWvTSWZ\nTW3cSiJWrxbp0aPdPnzTTd7t/didI+3ZAiKBgEhRkcjcuTKIJsfzjPF3XWOi7elq41bSDWrjVvKd\nVauiIyOPOMK7vW13Xr7cfRZsE7a8tLZCczPs3cuQI53ztSaaIUfmBd+8uT1SM7bYg6Z1VbKJEb81\nmZJgwoQJsmbNmrT3q3QdIu3axcXx5hEvYqutz58PtdMDnsm0t4+axNFb/xS1oFhSomKr5A/GmDdF\nZIKftjrjVrJKrF37rruSE21wmQW7iXYwCDNmMPCFFTpDVroMKtxK1mhosLxFNmywnq9aBTffnJxo\nOxHrGSJtW3MgRN0R9QT+cwnDaqzwdiezR7rRwBwl06hwK1njxRfh0CHrcTAIf/5zGjqtrubwD+ZE\n7bJD2wtaD3Hvru/7dv9LB/X1yUd5KkrS+FnBBM4F3gc+An6SqH138ir55JNPZNKkSdLkEY337LPP\nSjAYlEcffVQmTZokjY2NruesW7dO+vTpI+vXr4/r3+nxI488IsFgUCorK8P9OfUxduxYMcZISUlJ\neL997OSTT5aamhppbGyUk08+WY455hgB5LHHHosag328pqZGmpqa5OGHH45qZ1934cKFEggE5Nhj\nj5XS0tLw9VavFikuFoFnBYwAUlxcLCNGjJBAICDFxcWycuVKqa6ulpKSEunZs6esWLFCevbsKSUl\nJTJ+/Pjw6/3lL38pwWBQ7gfpDbISZBJIE8gnIDUgY0FKQKCHQEBgrJSWPifGlAqMlVCoRubPb39N\nCxYsEECKiopk1KhRUlpaKitWrJA+ffrE/V25cqXU1NRIdXW11NTUyHPPPSd9+vSRsrIVAn0EVgpU\nC5QIFEsweIwYY6S4uFiOOeYYCQQCMmrUKBk/fnz4va+urpZevXqFjx9//PHhY5GfF6fP0/jx4+PG\nW1paKmPHjpVRo0aF+xs/fryMHTtWevXqJePHjw/3bf9v7bZDhgwRQI455pjw6+vZs6f06tVL1q9f\n7/i5j/zcrVu3Lur/Fvsa7LYrV66Uk08+Ofw+2J/tmpqa8PPuAkl4lfgR7SDwV+BoIASsB0Z7ndOd\nhLuurk4CgYDU1dW5tunbt68AEgqFJBAISEVFhes5FRUVAkhFRUVc/06PQ6GQbRkI9+fUB+0WhPD+\n2GP2efYWCoWixhB5vK6uLnxtu519PBAIuF5v9WqRoqK+UccjN/u9at9CcWOKfB6wz2t7XAcyx6Vv\na4vtv/01xY478nqxf2PH2f485HId7y32vY89Fvl5cfo8uY3Xz3Vj+0j0f6moqHD83Ed+7mL7i30N\n9vHY99H+bMd+prsDyQh3Qq8SY8xE4DYROaft+c0AIvLPbud0B6+S4uJiDhyIz/dcVFTE/v37ATBu\nMdMO5zj1pSiKReT3qquSbq+SI4G/RTzf1rYv9qKzjTFrjDFrdu7c6W+knZiNGzfy/e9/n5I2p+KS\nkhJqa2vZtGlTuM2zzz4bPu5EcXFx+Jx169ZRXl4edbxnz54UtVWpDQaDBNuKKQaDQdebQo8ePTyf\n2/Tv35/JkyeH+/SisLAwYZtEDBkyhIULF7qOR0lMKBTy9f/KJiUlJXznO99hyJAhvtr36NHD8zsR\nSzAY5IILLoj6XilpXJwUkUUiMkFEJgwYMCBd3eYtgwcPpnfv3hw4cCA8Y+7duzeDIpLzT5kyxVWo\ngsEgBw8eDJ9TVVVFz549o9qEQiEOHTpEUVERLS0ttLS0hB+7/VKK/VK4fUkGDhzIcccdR0tsSfQY\nAoEALS0tHRaMPn368IMf/CCpL63STjAY5PDhwxkRbr+/DGOxP/fDhg2jd+/eCdsHg0Gam5uTunm3\ntLQwcODAqO+V4k+4PwaOing+pG1ft2f79u3MmTOH1157jTlz5vDpp5/Gtdm3bx99+/aloqKCUCiE\nMYaKigpOP/30uHM+//xzKioqWLFiBRUVFezbty/c//Dhwxk+fHjUY7u/vn370rt3bwKBAPv27Yvr\nIxCw/s1FRUXhL/7u3bvZvn07w4cPZ9q0aXzzm9+koMAqiGSMobS0FGMMra2tzJkzh9NPP53S0lKO\nPPJIpk2bFm533XXXEQqFAKioqAj3YR+3+9rdViZ9X2xKPVd6JnhuYX+A+wL9geFAiTHhXyrOhIDe\nbWcXUFLS/gMyVsTs12C/16FQiIKCAgoLCwkEAowYMYJAIEBhYSEFBQXh4+39JBbaHj16YIwJ/59i\nKS0tDX9eDh8+HPd5CgaDlJSUhJ9Hbm6vy/4sFBQUUFZWRmlpacJxFhUVhWfWJSUlUZ/7yM+u/Rkr\nKioK36gLCgrCr8H+jJ544okUFhZijGHEiBEceeSRlJSUUFJSwrRp05g2bRrDhw93/F51exIZwbHK\nm23E+k7Yi5MVXud0p8VJJb0kygFijMjXFLknLQGRmhrX/uvqrD7c+ncaj1OuknS9HkWxIZ1eJVZ/\nTAU+wPIumZeovQq30hG8xDIymZTjNniwZ6EEt0RV5eXO4ygsjG5XWJiaeHdE/JXuQTLCrblKlLwm\nNi/J1Klw34Lon/32J9h897tWHLuHPTTgktLEGCuiMpL+/WHXrvi2ZWXw2WdJvhBFSYDmKlG6BE5R\niD9YWB3XLlwE+KWXPEUbksuj7STaXvsVJVuocCt5S2x5sLVUM04anRtXVFjCnQCv4giK0llQ4Vby\nltgyYNW4iPa0afDcc1BZGXcoNuET+M+jXVbmfDm3/YqSLVS4lbzFdxmwxx6DSy+N2+2W8AmcswTG\nivy0aVb9y0hCIbjnntRej6KkCxVuJW/xU2sSsOzaK1bE7U6mEruTyC9dClddFT07f+ghzeGt5B4V\nbiVviS0PFkcwCDNnwtq1jouSyVRidxP5Rx9NetiKknFUuJW8JrLaTRhbsLdtg8WLXT1J3EwtgUB8\nfmw3kd+1S3NrK/mHCrfSeaiqgrlzEwq2jZuppaUlXoD92tPdTC2Kkk00AEfp0tTXw4wZlljHUl5u\nzebtdrNnx5tLnHAK1lGUjqIBOIrSRm2tu8hGmkdi7enl5e5uf769XRQlQ6hwK10ev9GSsdXj77lH\ng3WU/ESFW+nypBot6TQLdwvWUZRsUpC4iaJ0bmyhjUxWNX++PwGurVWhVvIPFW6lW6ACrHQl1FSi\nKIrSyVDhVhRF6WSocCuKonQyVLgVRVE6GSrciqIonYy0e5UYY2YDnxljtqS77zTQH8jXaoH5PDbI\n7/Hp2FJDx5YamRpbud+Gac9VYoxZ4zfePtvo2FInn8enY0sNHVtq5MPY1FSiKIrSyVDhVhRF6WRk\nQrgXZaDPdKFjS518Hp+OLTV0bKmR87FlJB+3oiiKkjnUVKIoitLJyIhwG2P+xRjznjHmLWPMU8aY\nIzJxnVQwxlxsjHnHGNNqjMmLVWtjzLnGmPeNMR8ZY36S6/FEYox5yBizwxjzdq7HEokx5ihjzAvG\nmA1t/8/rcj2mSIwxRcaY140x69vGd3uuxxSJMSZojFlnjPl9rscSizFmszHmL8aYRmNMXpXSMsYc\nYYx5vE3f3jXGTMzFODI1434OGCMilcAHwM0Zuk4qvA1cALyU64GA9QUCfgN8BxgNXGaMGZ3bUUWx\nBDg314Nw4DDwf0RkNFAD/DDP3reDwGQRGQdUAecaY2pyPKZIrgPezfUgPPi2iFTl2u3OgXuAP4rI\nKGAcOXoPMyLcIvKsiBxue/oaMCQT10kFEXlXRN7P9TgiOAn4SEQ2isgh4BHgH3I8pjAi8hKwO9fj\niEVEmkRkbdvjr7C+QEfmdlTtiMXetqeFbVteLCgZY4YAfwc8kOuxdCaMMX2AScCDACJySET25GIs\n2bBxzwL+Xxau01k5EvhbxPNt5JEAdQaMMcOAauDPuR1JNG3miEZgB/CciOTL+H4F3ATka8ljAVYa\nY95si8TOF4YDO4HFbWamB4wxPXMxkJSF2xiz0hjztsP2DxFt5mH9pK1Px2DTOTala2CM6QU8AVwv\nIl/mejyRiEiLiFRh/eI8yRgzJtdjMsacB+wQkTdzPRYPvtX2vn0HywQ2KdcDaqMAGA8sEJFq4Gsg\nJ2tSKecqEZGzvI4bY2YC5wFnSpZ9DhONLc/4GDgq4vmQtn1KAowxhViiXS8iT+Z6PG6IyB5jzAtY\nawW5XuQ9Ffh7Y8xUoAjobYxZLiKX53hcYUTk47a/O4wxT2GZE/NhTWobsC3il9Pj5Ei4M+VVci7W\nT7G/F5F9mbhGF+IN4FhjzHBjTAi4FPhdjseU9xhjDJat8V0R+bdcjycWY8wA25vKGFMMTAHey+2o\nQERuFpEhIjIM67P2fD6JtjGmpzGm1H4MnE3ub3YAiMinwN+MMSPbdp0JbMjFWDJl4/41UAo81+bS\nszBD10kaY8z5xphtwETgD8aY/8nleNoWcX8E/A/WAtujIvJOLscUiTHmYaABGGmM2WaMuSrXY2rj\nVGA6MLntM9bYNovMFwYDLxhj3sK6OT8nInnnepeHDAReMcasB14H/iAif8zxmCK5Fqhv+79WAXfl\nYhAaOakoitLJ0MhJRVGUToYKt6IoSidDhVtRFKWTocKtKIrSyVDhVhRF6WSocCuKonQyVLgVRVE6\nGSrciqIonYz/D5Bm+jZoIVUIAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x220eefe3128>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#!/usr/bin/env python\n",
    "# coding: utf-8\n",
    "#copyRight by heibanke\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def pca(dataMat, topNfeat=5):  \n",
    "    meanVals = np.mean(dataMat, axis=0)  \n",
    "    meanRemoved = dataMat - meanVals                     #减去均值  \n",
    "    covMat = np.cov(meanRemoved, rowvar=0)               #求协方差方阵  \n",
    "    eigVals, eigVects = np.linalg.eig(covMat)            #求特征值和特征向量  \n",
    "    eigValInd = np.argsort(eigVals)[:-(topNfeat + 1):-1] #对特征值进行排序  \n",
    "    redEigVects = eigVects[:, eigValInd]                 #除去不需要的特征向量  \n",
    "    lowDDataMat = meanRemoved.dot(redEigVects)           #求新的数据矩阵  \n",
    "    reconMat = (lowDDataMat.dot(redEigVects.T)) + meanVals\n",
    "    reduceData = lowDDataMat + np.mean(dataMat)\n",
    "\n",
    "    return reduceData,reconMat\n",
    "    \n",
    "def optimal(A,b):\n",
    "    B = A.T.dot(b)\n",
    "    AA = np.linalg.inv(A.T.dot(A))#求A.T.dot(A)的逆\n",
    "    P=AA.dot(B)\n",
    "    return A.dot(P)    \n",
    "    \n",
    "N=200\n",
    "x=np.linspace(2,4,N)\n",
    "y=x*2-4\n",
    "\n",
    "x1=x+(np.random.rand(N)-0.5)*1.5\n",
    "y1=y+(np.random.rand(N)-0.5)*1.5\n",
    "\n",
    "data = np.array([x1,y1])\n",
    "a,b=pca(data.T,1)\n",
    "\n",
    "one=np.ones((len(x1),1))#len(x)得到数据量\n",
    "xx=x1.reshape((len(x1),1))\n",
    "A=np.hstack((xx, one))#两个100x1列向量合并成100x2,(100, 1) (100,1 ) (100, 2)\n",
    "yy=optimal(A,y1)\n",
    "plt.plot(xx,yy,color='b',linestyle='-',marker='.',label=u\"nihe\")\n",
    "\n",
    "plt.plot(x,y,color='g',linestyle='-',marker='',label='ideal') \n",
    "plt.plot(x1,y1,color='b',linestyle='',marker='o',label='noise')\n",
    "plt.plot(b[:,0],b[:,1],color='r',linestyle='',marker='>',label='recon')\n",
    "plt.plot(a[:,0],np.zeros(N),color='k',linestyle='',marker='*',label='lowD')\n",
    "\n",
    "plt.legend()\n",
    "plt.axis('equal')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[3 1]\n",
      " [2 4]]\n",
      "[ 2.  5.]\n",
      "[[-0.70710678 -0.4472136 ]\n",
      " [ 0.70710678 -0.89442719]]\n",
      "[[-0.52573111 -0.85065081]\n",
      " [-0.85065081  0.52573111]]\n",
      "[ 5.11667274  1.95439508]\n",
      "[[-0.64074744 -0.76775173]\n",
      " [-0.76775173  0.64074744]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "A=np.array([[3,1],[2,4]])\n",
    "eigVals, eigVects = np.linalg.eig(A)\n",
    "print(A)\n",
    "print(eigVals)\n",
    "print(eigVects)\n",
    "u,s,v = np.linalg.svd(A)\n",
    "print(u)\n",
    "print(s)\n",
    "print(v)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
